{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load the required packages\n",
    "using ODE\n",
    "using JLD\n",
    "using ForwardDiff\n",
    "\n",
    "# Set precision to quadruple precision (QP)\n",
    "# QP has 113 bits of precision\n",
    "set_bigfloat_precision(113);\n",
    "\n",
    "# Number of subdivisions of the rope\n",
    "global n = 40;\n",
    "\n",
    "# Define the right-hand function\n",
    "function rope(t,x)\n",
    "    \n",
    "    n2 = n*n; # n^2\n",
    "    n3by4 = convert(Int64,3*n/4); # 3*n/4\n",
    "    \n",
    "    # Force in x-direction\n",
    "    Fx = parse(BigFloat,\"0.4\");\n",
    "    # Force in y-direction\n",
    "    Fy = cosh(4*t-2.5)^(-4);\n",
    "    \n",
    "    # Compute required matrices\n",
    "    c = -cos(x[1:n-1]-x[2:n]);\n",
    "    cDiag = [one(BigFloat);\n",
    "        2*ones(BigFloat,n-2);\n",
    "        3*one(BigFloat)];\n",
    "    C = spdiagm((c,cDiag,c),(-1,0,1));\n",
    "    \n",
    "    d = -sin(x[1:n-1]-x[2:n]);\n",
    "    D = spdiagm((-d,d),(-1,1));\n",
    "    \n",
    "    # Compute the inhomogeneous term\n",
    "    v = -(n2+n/2-n*[1:n;]).*sin(x[1:n])-n2*sin(x[1:n])*Fx;\n",
    "    v[1:n3by4] = v[1:n3by4] + n2*cos(x[1:n3by4])*Fy; \n",
    "    \n",
    "    w = D*v+x[n+1:2*n].^2;\n",
    "    u = C\\w;\n",
    "    \n",
    "    # Write down the system\n",
    "    return [x[n+1:2*n],C*v + D*u];\n",
    "end\n",
    "\n",
    "# Set up the initial conditions\n",
    "t0 = zero(BigFloat);\n",
    "T = parse(BigFloat,\"3.723\");\n",
    "x0 = zeros(BigFloat,2*n);\n",
    "\n",
    "# Set the tolerance \n",
    "Tol = parse(BigFloat,\"1e-33\");\n",
    "\n",
    "# Solve and get the solution at T = tEnd\n",
    "(t_rope,x_tmp_rope) = ode78(rope,x0,[t0;T];\n",
    "reltol=Tol,abstol=Tol,points=:specified);\n",
    "\n",
    "x_ref = Array{BigFloat}(n);\n",
    "\n",
    "x_ref[:] = x_tmp_rope[2,1][1:n];\n",
    "\n",
    "save(\"refSolRope.jld\",\"x_ref\",x_ref);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
